## Created by EmpowerStats @ Thu, 05 Dec 24 19:29:57 +0800## 
#******************** Regarding ALL Following R Code   ******************** 
#*** COPYRIGHT (c) 2010, 2021 X&Y Solutions, ALL RIGHT RESERVED *********** 
#******************* www.EmpowerStats.com ********************************* 
#************************************************************************** 
Sys.setlocale(category = 'LC_ALL', locale = 'English_United States.1252'); 
.libPaths(file.path(R.home(),'library')); 
library(doBy); 
options(timeout=600); 
library(plotrix); 
library(stringi); 
library(stringr); 
library(survival); 
library(rms); 
library(nnet); 
library(car); 
library(mgcv); 
pdfwd<-6; pdfht<-6; 
load('/Users/liukangning/Desktop/EmpowerROS/Analysis/cognitivedepression/datamergeall.Rdata'); 
if (length(which(ls()=='EmpowerStatsR'))==0) EmpowerStatsR<-get(ls()[1]); 
names(EmpowerStatsR)<-toupper(names(EmpowerStatsR)); 
originalVNAME<-names(EmpowerStatsR); 
ofname<-'cognitivedepression_7_tbl'; 
vname<-c(NA,'EDUCATION','EDUCATION.1','EDUCATION.2','EDUCATION.3','EDUCATION.4','EDUCATION.5','EDUCATION.9','MARITAL','MARITAL.1','MARITAL.2','MARITAL.3','MARITAL.4','MARITAL.5','MARITAL.6','MARITAL.77','BP','BP.1','BP.2','BP.9','CORONARY','CORONARY.1','CORONARY.2','CORONARY.9','STROKE','STROKE.1','STROKE.2','STROKE.9')[-1]; 
vlabel<-c(NA,'受教育程度','  1','  2','  3','  4','  5','  9','婚姻状况','  Married','  Widowed','  Divorced','  Separated','  Never married','  Living with partner','  Refused','高血压','  1','  2','  9','冠心病','  1','  2','  9','中风','  1','  2','  9')[-1]; 
varused4this <- c('EDUCATION','MARITAL','BP','CORONARY','STROKE'); 
pkgs<-c('MASS','Matrix','nnet','car','arm','lattice','lme4','R2WinBUGS','coda','abind','mi','missForest'); 
for (g in pkgs) {  
if (!(g %in% rownames(installed.packages()))) install.packages(g,repos='https://cloud.r-project.org'); 
}
library(MASS); 
library(Matrix); 
library(nnet); 
library(car); 
library(arm); 
library(lattice); 
library(lme4); 
library(R2WinBUGS); 
library(coda); 
library(abind); 
library(mi); 
library(missForest); 
WD <- EmpowerStatsR; rm(EmpowerStatsR); gc(); 
title<-'缺失数据多重插补 (MI)'; 
wd.subset=''; 
xvname<-c('MARITAL','EDUCATION','CORONARY','STROKE','BP'); 
sxf<-c(0,0,0,0,0); 
xlv<-c(7,6,3,3,3); 
par2<-5; 
chk<- 0; 
dec<-4; 
##R package##  MASS Matrix nnet car arm lattice lme4 R2WinBUGS coda abind mi missForest ##R package##;
printTxt<-function(txt,cname,rname) {
  tmp<-as.matrix(txt,ncol=1); rname1<-rname
  if (length(rname1)==1) rname1=rep(rname,times=nrow(tmp))
  colnames(tmp)=cname;  rownames(tmp)=rname1;print(tmp,quote=F)
}
mat2htmltable<-function(mat, addrcname=FALSE) {
  if (addrcname) {mat<-cbind(c(" ",rownames(mat)),rbind(colnames(mat),mat)); }
  t1<- apply(mat,1,function(z) paste(z,collapse="</td><td>"))
  t2<- paste("<tr><td>",t1,"</td></tr>")
  return(paste(t2,collapse=" "))
}
tmp.nimp<-par2
if (!is.numeric(tmp.nimp) | tmp.nimp<1) tmp.nimp<-5
if (tmp.nimp>30) tmp.nimp<-30
nx<-length(xvname)
tmp.d<-WD[,xvname]
tmp.mdf<-missing_data.frame(tmp.d)
png(paste(ofname,"_missing.png",sep=""),width=720,height=720);try(image(tmp.mdf));dev.off()
png(paste(ofname,"_hist.png",sep=""),width=720,height=720);hist(tmp.mdf);dev.off()
v.type<-c("binary","ordered-categorical","unordered-categorical","count","proportion","continuous","bounded-continuous","positive-continuous","nonnegative-continuous","semi-continuous");
for (i in (1:nx)) {
  if (sxf[i]>=1 & sxf[i]<=length(v.type)) {
     tmp.mdf<-change(tmp.mdf,y=xvname[i],what="type",to=v.type[sxf[i]])
  }
}
tmp.oth<-NA
if (length(xvname)<length(colnames(WD))) tmp.oth<-WD[, -match(xvname,colnames(WD))]
tmp.dd<-cbind(mi.iter=0,WD[,xvname],tmp.oth)
okmissForest<-1;
if (chk==1) {
  for (v in xvname[xlv>0]) tmp.d[,v]<-factor(tmp.d[,v])
  mistmp<-try(missForest(tmp.d, maxiter=10, ntree=50, parallelize='no'));
  if (substr(mistmp[1],1,5) != "Error") {
	tmp.dd<-rbind(tmp.dd, cbind(mi.iter=1, mistmp$ximp[,xvname], tmp.oth))
    for (i in 2:tmp.nimp) {
      mistmp<-missForest(tmp.d, maxiter=10, ntree=50, parallelize='no')
      tmp.dd<-rbind(tmp.dd,cbind(mi.iter=i,mistmp$ximp[,xvname],tmp.oth))
    }
  } else {okmissForest<-0;}
} else {okmissForest<-0;}
if (okmissForest==0) {
  tmp.imp<-mi(tmp.mdf,n.iter=30,n.chains=4,max.minutes=20,seed=NA,verbose=TRUE,parallel=.Platform$OS.type!="windows")
  tmp.cc<-complete(tmp.imp,30)
  tmp.dd<-rbind(tmp.dd, cbind(mi.iter=1,tmp.cc[[1]][,xvname], tmp.oth))
  if (tmp.nimp>1) {for (i in 2:tmp.nimp) {tmp.dd<-rbind(tmp.dd,cbind(mi.iter=i,tmp.cc[[i]][,xvname], tmp.oth))}}
}
write.table(tmp.dd,file=paste(ofname,"_mi.xls",sep=""),quote=FALSE,col.names=TRUE,row.names=FALSE,sep="\t")
tb<-table(tmp.mdf@patterns)
tb<-rbind(c("Pattern","N"),cbind(names(tb),tb))
ts0<-round(t(apply(WD[,xvname],2,function(x) summary(as.numeric(x[!is.na(x)])))),2)
ts1<-round(t(apply(tmp.dd[tmp.dd[,1]>0,xvname],2,function(x) summary(as.numeric(x)))),2)
ts0<-cbind(c("",rownames(ts0)),rbind(colnames(ts0),ts0))
ts1<-cbind(c("",rownames(ts1)),rbind(colnames(ts1),ts1))
w<-c("<!DOCTYPE html><html lang='zh'><head><meta charset='utf-8'></head><body>")
w<-c(w,"<h2>Multiple imputation</h2>")
w<-c(w,"Missing data patterns:</br><table border=3>", mat2htmltable(tb), "</table>")
w<-c(w,"</br>Original data summary:</br><table border=3>", mat2htmltable(ts0), "</table>")
w<-c(w,"</br>Imputed data summary:</br><table border=3>", mat2htmltable(ts1), "</table>")
if (okmissForest==0) {
  w<-c(w,"Using R mi package");
  w<-c(w,"</br>If imputed data range out of limits, try change variable type and re-do imputation.")
  w<-c(w,"</br>Variable type options:")
  w<-c(w,paste("</br>",(1:length(v.type)),"=",v.type))
} else {
  w<-c(w,"Using random forest imputation (missForest package)");
  w<-c(w,paste("</br></br>Created by EmpowerStats (www.empowerstats.com) and R on",Sys.Date()))
  w<-c(w,"</body></html>")
}
fileConn<-file(paste(ofname,".htm",sep="")); writeLines(w, fileConn)

if (okmissForest==0) {
  sink(paste(ofname,".htm",sep=""), append=TRUE)
  nr<-nx*2+10
  printTxt(paste("</br><textarea rows=", nr, " cols=100>",sep=""),"","")
  show(tmp.mdf);
  printTxt("</textarea></br>","","")
  printTxt(wd.subset,"","")
  printTxt(paste("</br>Created by EmpowerStats (www.empowerstats.com) and R on",Sys.Date()),"","")
  printTxt("</body></html>","","")
  sink()
}